i = 0
nn.results <- list()
size_seq <- seq(4,12,2)
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Fit NN
  nn.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("nnet",
                                       "pROC")) %dopar% {   # Loop over CV runs
                 set.seed(52*cv) 
                 shuffle <- sample(1:N,
                                   N,
                                   replace=FALSE)
                 pred.prob <- matrix(NA,nrow=N,ncol=length(size_seq))
                 sampleSplit = c(0,round((N/5*c(1:5))))
                 for (f in 1:5) {
                   # Remove constant variables from predictors
                   train.RHS <- dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                 rhs]
                   test.RHS <- dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                              rhs]
                   not.constant <- apply(train.RHS,2,sd)!=0
                   train.RHS <- train.RHS[,not.constant]
                   test.RHS <- test.RHS[,not.constant]
                   # Find PC weights
                   pcs.fit <- prcomp(train.RHS,center = TRUE, scale = TRUE)
                   train.pcs <- predict(pcs.fit,
                                        train.RHS)[,1:min(nn_params$pcNum[[country]],
                                                          ncol(train.RHS))]
                   test.pcs <- predict(pcs.fit,
                                       test.RHS)[,1:min(nn_params$pcNum[[country]],
                                                        ncol(train.RHS))]
                   for (s in 1:length(size_seq)) {
                     # Fit Neural Network
                     cv.fit = nnet(x = train.pcs,
                                   y = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                               dv]),
                                   size = size_seq[s], 
                                   decay = nn_params$decay,
                                   trace = nn_params$trace
                     )
                     pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                               s] <- predict(cv.fit,
                                            newdata = as.matrix(test.pcs))
                   }
                 }
                 cv_results <- c("CV Run"=0,
                                 "Size"=0,
                                 "MSE"=100)
                 # Get error rates for each lambda used
                 for (s in 1:length(size_seq)) {
                     mse <- mean((dta.past[,dv] - pred.prob[,s])^2)
                     if (mse<cv_results["MSE"]) {
                       cv_results <- c("CV Run"=cv,
                                       "Size" = size_seq[s],
                                       "MSE" = mse)
                       best.s <- s
                     }
                 }
                 cv_results<- c(cv_results,
                                "DT" = as.numeric(mostaccDT(y=dta.past[,dv],
                                                            yhat=pred.prob[,s],
                                                            J=1000)),
                                "DTsens" = as.numeric(bestDT(y=dta.past[,dv],
                                                            yhat=pred.prob[,s],
                                                            J=1000))
)
                 return(cv_results)
  }
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  nn.results[[i]]$param_dist <- cv_results
  if (cv.runs > 1) {
    nn.results[[i]]$size <- round(median(cv_results[,"Size"]))
    nn.results[[i]]$dt <- mean(cv_results[,"DT"])
    nn.results[[i]]$dtsens <- mean(cv_results[,"DTsens"])
  } else {
    nn.results[[i]]$size <- cv_results["Size"]
    nn.results[[i]]$dt <- cv_results["DT"]
    nn.results[[i]]$dtsens <- cv_results["DTsens"]
  }
 
  dta.past.rhs <- dta.past[,rhs]
  not.constant <- apply(dta.past.rhs,2,sd)!=0
  dta.past.rhs <- dta.past.rhs[,not.constant]
  pcs.fit <- prcomp(dta.past.rhs,
                    center = TRUE, 
                    scale = TRUE)
  insample.pcs <- predict(pcs.fit,dta.past.rhs)[,1:min(nn_params$pcNum[[country]],
                                                       ncol(dta.past.rhs))]
  outsample.pcs <- predict(pcs.fit,
                           dta[dta[,t.var]== year,rhs])[,1:min(nn_params$pcNum[[country]],
                                                               ncol(dta.past.rhs))]
  set.seed(i)
  mod <- nnet(x = insample.pcs,
                   y = as.matrix(dta[dta[,t.var]< year,
                                          dv]),
                   size = nn.results[[i]]$size,
                   decay = nn_params$decay,
                   trace = nn_params$trace
  )
  nn.results[[i]]$fit.oos <- predict(mod,
                                      newdata = outsample.pcs)
  nn.results[[i]]$actual.oos <- as.matrix(dta[dta[,t.var] == year,
                                                    dv])
  print(year)
}
save(nn.results,file=paste(modeldir,"/",
                              country,
                              "_nn_",
                              v,
                              "_",
                              rhs.group,
                              "_mse",
                              ".RData",
                              sep = ""))